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Abstract The escape dynamics in a simple analytical gravi¬ 
tational model which describes the motion of stars in a Seyfert 
galaxy is investigated in detail. We conduct a thorough nu¬ 
merical analysis distinguishing between regular and chaotic 
orbits as well as between trapped and escaping orbits, con¬ 
sidering only unbounded motion for several energy levels. 
In order to distinguish safely and with certainty between or¬ 
dered and chaotic motion, we apply the Smaller ALingment 
Index (SALI) method. It is of particular interest to locate 
the escape basins through the openings around the collinear 
Lagrangian points Li and L 2 and relate them with the cor¬ 
responding spatial distribution of the escape times of the or¬ 
bits. Our exploration takes place both in the physical (v,y) 
and in the phase (v,i) space in order to elucidate the es¬ 
cape process as well as the overall orbital properties of the 
galactic system. Our numerical analysis reveals the strong 
dependence of the properties of the considered escape basins 
with the total orbital energy, with a remarkable presence of 
fractal basin boundaries along all the escape regimes. It was 
also observed, that for energy levels close to the critical es¬ 
cape energy the escape rates of orbits are large, while for 
much higher values of energy most of the orbits have low 
escape periods or they escape immediately to infinity. We 
also present evidence obtained through numerical simula¬ 
tions that our model can describe the formation and the evo¬ 
lution of the observed spiral structure in Seyfert galaxies. We 
hope our outcomes to be useful for a further understanding 
of the escape mechanism in galaxies with active nuclei. 
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1 Introduction 

Over the last decades a huge amount of research work has 
been devoted on the subject of escaping particles from open 
dynamical systems. Especially the issue of escape in Hamil¬ 
tonian systems is a classical problem in nonlinear dynam- 
ics (e.g., (28l|32|-|35l[88l). It is well known, that several 
types of Hamiltonian systems have a finite energy of es¬ 
cape. For values of energy lower than the escape energy the 
equipotential surfaces of the systems are close which means 
that orbits are bound and therefore escape is impossible. For 
energy levels above the escape energy on the other hand, 
the equipotential surfaces open and exit channels emerge 
through which the particles can escape to infinity. The lit¬ 
erature is replete with studies of such “open" Hamiltonian 
systems (e.g., tSl 1^ [78l [88l IMl - ITOOl [1^91 IT^ At this 
point we should emphasize, that all the above-mentioned 
references on escapes in Hamiltonian system are exemplary 
rather than exhaustive, taking into account that a vast quan¬ 
tity of related literature exists. 

Nevertheless, the issue of escaping orbits in Hamilto¬ 
nian systems is by far less explored than the closely related 
problem of chaotic scattering. In this situation, a test particle 
coming from infinity approaches and then scatters off a com¬ 
plex potential. This phenomenon is well investigated as well 
interpreted from the viewpoint of chaos theory (e.g., M- 

-[ 95 I). Chaotic scattering has also been applied in the astro- 
physical context of many aspects such as scattering off black 
holes (e.g., EIEtI) and three-body stellar systems (e.g., Col 
cniioiisii). The related invariant manifolds of the chaotic 
saddle is directly associated with the chaotic dynamical be¬ 
havior (e.g., ISTJ). In particular, the chaotic saddle is defined 
as the intersection of its stable and unstable manifolds Ell, 
while hyperbolic and non-hyperbolic chaotic saddles may 
occur in dynamical systems (e.g., Col). More details on the 
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issue of chaotic scattering and escape from chaotic systems 
can be found to the recent reviews ll96ll and ii , respectively. 

In open Hamiltonian systems an issue of paramount im¬ 
portance is the determination of the basins of escape, simi¬ 
lar to basins of attraction in dissipative systems or even the 
Newton-Raphson fractal structures. An escape basin is de¬ 
fined as a local set of initial conditions of orbits for which 
the test particles escape through a certain exit in the equipo- 
tential surface for energies above the escape value. Basins 
of escape have been studied in many earlier papers (e.g., 
113 EH EH [13). The reader can find more details regard¬ 
ing basins of escape in EH, while the review EH provides 
information about the escape properties of orbits in a multi¬ 
channel dynamical system of a two-dimensional perturbed 
harmonic oscillator. The boundaries of an escape basins may 
be fractal (e.g., El [13) or even respect the more restrictive 
Wada property (e.g., IH), in the case where three or more 
escape channels coexist in the equipotential surface. 

One of the most characteristic Hamiltonian systems of 
two degrees of freedom with escape channels is undoubtedly 
the well-known Henon-Heiles system (491 . A huge load of 
research on the escape properties of this system has been 
conducted over the years (e.g., |[I]-[3l[7l[8l). Escaping orbits 
in the classical Restricted-Three-Body-Problem (RTBP) is 
another typical example (e.g., 17^1801110411 ). Furthermore, 
escaping and trapped motion of stars in stellar systems is an 
another issue of great importance. In a recent article 11071 . 
we explored the nature of orbits of stars in a galactic-type 
potential, which can be considered to describe local motion 
in the meridional (R,z) plane near the central parts of an 
axially symmetric galaxy. It was observed, that apart from 
the trapped orbits there are two types of escaping orbits, 
those which escape fast and those which need to spend vast 
time intervals inside the equipotential surface before they 
find the exit and eventually escape. Furthermore, the chaotic 
dynamics within a star cluster embedded in the tidal field 
of a galaxy was explored in l43l . In particular, by scanning 
thoroughly the phase space and obtaining the basins of es¬ 
cape with the respective escape times it was revealed, that 
the higher escape times correspond to initial conditions of 
orbits near the fractal basin boundaries. 

The aim of this work is to numerically investigate the 
escape dynamics in a Seyfert galaxy with a rotating active 
central nucleus. The recent paper of Ernst & Peters (2014) 
m is used as a guide since we follow the same numeri¬ 
cal methods and computational approach. The present paper 
however, contains several new aspects with respect to El 
that can be considered as further developments in the field of 
open galactic systems. In particular: (i) we proceed one step 
further by using the SAFI method in order to distinguish be¬ 
tween regular and chaotic trapped motion, (ii) we display 
the evolution of the percentages of all types of orbits in both 
the physical and phase space as a function of the value of the 


energy integral, (iii) we conduct an overview analysis thus 
revealing how the orbital structure of the dynamical system 
is infiuenced by the total energy when we use a continu¬ 
ous spectrum of energy values rather than some few discrete 
levels and (iv) we relate the escape times of orbits with the 
energy as well as with the width of the symmetrical escape 
channels. 

The article is organized as follows: In Sectionj^we present 
in detail the structure and the properties of our Seyfert galaxy 
model. All the computational methods we used in order to 
determine the character of orbits are described in Section O 
In the following Section, we conduct a thorough numerical 
investigation revealing the overall orbital structure (bounded 
regions and basins of escape) of the galaxy and how it is af¬ 
fected by the total orbital energy. In Section we show that 
our dynamical model can simulate, in away, the creation as 
well as the evolution of spiral arms of a real Seyfert galaxy. 
Our paper ends with Sectionj^ where the discussion and the 
conclusions of this research are given. 


2 Theory of the galaxy model 


For the description of the Seyfert galaxy we use the follow¬ 
ing mass-model potential 




-GM, 

■s/x^ -\- by^ + ' 


( 1 ) 


where G is the gravitational constant. Mg is the total mass of 
the galaxy, Z? is a parameter leading to deviation from axial 
and spherical symmetry, while c is the scale length of the 
galaxy which acts also as a softening parameter. The gener¬ 
alized anharmonic Plummer potential Q has been success¬ 
fully applied several times in the past in order to realisti¬ 
cally model non axially symmetric galaxies (see e.g., EH 
[73ll^rT08l ). At this point we would like to stress out that 
the single-component potential of Eq. ^ models the prop¬ 
erties of the Seyfert galaxy as a whole rather than an isolate 
component (e.g., the galactic bar). 

We consider the case where the dense and compact ac¬ 
tive nucleus of the Seyfert galaxy rotates clockwise at a con¬ 
stant angular velocity ^2 > 0. Therefore, the effective poten¬ 
tial is given by 


^ef{{x,y) = ^g{x,y) - 


~Y 


(x^+y^), 


( 2 ) 


which is in fact a cut through the equatorial plane of the 
three-dimensional (3D) potential. 

In our study, we use a system of galactic units, where the 
unit of length is 1 kpc, the unit of mass is 2.325 x IO^Mq 
and the unit of time is 0.9778 x 10^ yr (approximately 100 
Myr). The velocity unit is 10 km/s, while G is equal to unity 
(G = 1). The energy unit (per unit mass) is 100 km^s“^. 
In these units, the values of the involved parameters are: 
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Mg = 5000 (corresponding to 1.1625 xlO^^ M©), b = 2.5, 
c = 1.2, ^2 = 1.25 and they remain constant throughout the 
numerical analysis. The values of Mg and c were chosen hav¬ 
ing in mind the Seyfert galaxy NGC 1566 (see at the end of 
section!^ for more details). 

This galactic model has five equilibria called Lagrangian 
points II103II at which 

ox oy 

The isolines contours of constant effective potential as well 
as the position of the five Lagrangian points L/, / = 1,5 are 
shown in Fig.[^. Three of them, Li, L 2 , and L 3 , known as 
the collinear points, are located in the v-axis. The central 
stationary point L 3 at (x,y) = ( 0 , 0 ) is a local minimum of 
At the other four Lagrangian points it is possible for 
the test particle to move in a circular orbit, while appearing 
to be stationary in the rotating frame. For this circular or¬ 
bit, the centrifugal and the gravitational forces precisely bal¬ 
ance. The stationary points Li and L 2 at (v,y) = (±rL,0) = 
(±14.687185207778671,0) are saddle points, where rL is 
called Lagrangian radius. Let Li located at v = — while 
±2 be at V = ±rz^. The points L 4 and L 5 on the other hand, 
which are located at (x,y) = (0, ±12.626321712993062) are 
local maxima of the effective potential, enclosed by the banana¬ 
shaped isolines. The annulus bounded by the circles through 
Li, ±2 and L 4 , L 5 is known as the “region of coroation" (see 
also d). 

The equations of motion in the rotating frame are de¬ 
scribed in vectorial form by 

f =-V4>eff-2(i2 xr), (4) 


where r = (v,y,z) is the position vector, ^2 = (0,0, ^2) is the 
constant rotation velocity vector around the vertical z-axis, 
while the term —2 {Q x f) represents the Coriolis force. De¬ 
composing Eq. 0 into rectangular Cartesian coordinates 
(v,y), we obtain 




'eff 


= 


dx 

dy 


2ny, 


— 2^2i, 


(5) 


where the dot indicates derivative with respect to the time. 

In the same vein, the equations describing the evolution 
of a deviation vector w = (Sx,Sy,Sx,Sy) which joins the 
corresponding phase space points of two initially nearby or¬ 
bits, needed for the calculation of standard chaos indicators 
(the S ALI in our case) are given by the following variational 
equations 
(Sx) = 5x, 

(Sy) = Sy, 


(Sx) 


d^<P, 


eff 


dx^ 


Sx- 


d^<P, 


eff 


(5y) = - 


d^(p, 


eff 


dydx 


Sx — 


dxdy 

dy^ 


5y±2X2 5y, 


5y —2X2 Sx. 


( 6 ) 


Consequently, the corresponding Hamiltonian (known 
also as the Jacobian) to the effective potential given in Eq. 

Q reads 

Hj(x,y,x,y) = 1 (x^+f) + j) = E, (7) 

where x and y are the momenta per unit mass, conjugate to v 
and y, respectively, while E is the numerical value of the Ja¬ 
cobian, which is conserved since it is an isolating integral of 
motion. Thus, an orbit with a given value for it’s energy in¬ 
tegral is restricted in its motion to regions in which E < ^eff, 
while all other regions are forbidden to the star. The numer¬ 
ical value of the effective potential at the two Lagrangian 
points Li and ± 2 , that is <J>eff(~^L, 0 ) and <J>eff(^ l, 0 ), respec¬ 
tively yields to a critical Jacobi constant Cl = —507.828303111454, 
which can be used to define a dimensionless energy param¬ 
eter as 


^_ Cl-E 
Cl ’ 


( 8 ) 


where E is some other value of the Jacobian. The dimension¬ 
less energy parameter C makes the reference to energy levels 
more convenient. For E = Cl the equipotential surface en¬ 
closes the critical are^ while for a Jacobian value E > Cl, 
or in other words when C > 0, the equipotential surface is 
open and consequently stars can escape from the galaxy. In 
Fig.[^ we present the contour of ^Qf[(x,y) = E when C = 
0.01. We observe, the two openings (exit channels) through 
which the particles can leak out. In fact, we may say that 
these two exits act as hoses connecting the interior region 
(green color) of the galaxy where —rL<x<rL with the 
“outside world" of the exterior region (amber color). Exit 
channel 1 at the negative v-direction corresponds to escape 
through Lagrangian point Li , while exit channel 2 at the pos¬ 
itive v-direction indicates escape through ± 2 . The forbidden 
regions of motion within the banana-shaped isolines around 
±4 and ±5 points are shown in Fig.[^ with gray. It is evident, 
that for E < Cl the two forbidden regions merge together 
thus escape is impossible since the interior region cannon 
communicate with the exterior area. 


3 Computational methods 

In order to explore the escape dynamics of stars in the Seyfert 
galaxy model, we need to define samples of initial condi¬ 
tions of orbits whose properties (bounded or escaping mo¬ 
tion) will be identified. The best method for this purpose, 
would have been to choose the sets of initial conditions of 
the orbits from a distribution function of the system. This 
however, is not available, so we define for each value of the 
total orbital energy (all tested energy levels are above the 

^ In three dimensions the last closed equipotential surface of the ef¬ 
fective potential passing through the Lagrangian points L\ and L 2 en¬ 
closes a critical volume. 
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Fig. 1 (a-left): The isolines contours of the constant effective potential and the position of the five Lagrangian points, (b-right): The contour of the 
effective potential <pQff{x^y) = E when C = 0.01. The escaping orbits leak out through exit channels 1 and 2 of the equipotential surface passing 
either L\ or L 2 , respectively. The interior galactic region is indicated in green, the exterior region is shown in amber color, while the forbidden 
regions of motion are marked with grey. The two unstable Lyapunov orbits around the Lagrangian points are visualized in red. 


escape energy), dense grids of initial conditions regularly 
distributed in the area allowed by the value of the energy. 
Our investigation takes place both in the physical (x, y) and 
in the phase (x,x) space for a better understanding of the 
escape process. In both cases, the step separation of the ini¬ 
tial conditions along the axes (or in other words, the den¬ 
sity of the grids) is controlled in such a way that always 
there are about 10^ orbits (a maximum grid of 318 x 318 
equally spaced initial conditions of orbits), depending on the 
value of the Jacobian integral. Furthermore, the grids of ini¬ 
tial conditions of orbits whose properties will be examined 
are defined as follows: For the physical (x,y) space we con¬ 
sider orbits with initial conditions (xo,yo) with yo = 0, while 
the initial value of xq is always obtained from the energy in¬ 
tegral 0 as x'o = x(xo, yo, To, > 0- Similarly, for the phase 
(x,x) space we consider orbits with initial conditions (xo,xo) 
with yo = 0, while this time the initial value of y'o is obtained 
from the Jacobi integral 0. 

The equations of motion as well as the variational equa¬ 
tions for the initial conditions of all orbits were integrated 
using a double precision Bulirsch-Stoer FORTRAN 77 algo¬ 
rithm (e.g., 1841 ) with a small time step of order of 10“^, 
which is sufficient enough for the desired accuracy of our 
computations. Here we should emphasize, that our previous 
numerical experience suggests that the Bulirsch-Stoer inte¬ 
grator is both faster and more accurate than a double preci¬ 
sion Runge-Kutta-Fehlberg algorithm of order 7 with Cash- 
Karp coefficients. Throughout all our computations, the Ja¬ 
cobian energy integral (Eq. (0) was conserved better than 


one part in 10 ^ ^ although for most orbits it was better than 
one part in 10“^^. 

In dynamical systems with escapes an issue of paramount 
importance is the determination of the position as well as 
the time at which an orbit escapes. For all values of en¬ 
ergy smaller than the critical value Cl (escape energy), the 
equipotential surface is closed and therefore all orbits are 
bounded. On the other hand, when E > Cl the equipoten¬ 
tial surface is open and extend to infinity. The value of the 
energy itself however, does not furnish a sufficient condi¬ 
tion for escape. An open equipotential curve consists of two 
branches forming channels through which an orbit can es¬ 
cape to infinity (see Fig. [^). It was proved 1^ that in 
Hamiltonian systems there is a highly unstable periodic or¬ 
bit at every opening close to the line of maximum potential 
which is called a Lyapunov orbit. In fact, there is a family of 
unstable Lyapunov orbits around each collinear Lagrangian 
point 1761 . Such an orbit reaches the Zero Velocity Curv^ 
(ZVC), on both sides of the opening and returns along the 
same path thus, connecting two opposite branches of the 
ZVC. Lyapunov orbits are very important for the escapes 
from the system, since if an orbit intersects any one of these 
orbits with velocity pointing outwards moves always out¬ 
wards and eventually escapes from the system without any 
further intersections with the surface of section (e.g., 1 ^ . 
When E = Cl the Lagrangian points exist precisely but when 
£■ > Cl an unstable Lyapunov periodic orbit is located close 


^ The boundaries of the accessible regions of motion are the called 
Zero Velocity Curves, since they are the locus in the (x,y) space where 
kinetic energy vanishes. 
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Fig. 2 An orbit with initial conditions outside the interior region which 
however does not escape immediately from the system. In particular, 
even though the orbit is initiated outside but very close to L\ it eventu¬ 
ally escapes from the opposite exit channel, through L 2 . 

to each of these two points (e.g., (471). The two unstable 
Lyapunov orbits around Li and L 2 are shown in Fig. in 
red. 

These two Lagrangian points are saddle points of the ef¬ 
fective potential, so when £” > Cl, a star must pass close 
enough to one of these points in order to escape. Thus, in 
our galactic system the escape criterion is purely geometric. 
In particular, escapers are defined to be those stars moving in 
orbits beyond the Lagrangian radius {ri) thus passing one of 
the two Lagrangian points {L\ or L 2 ) with velocity pointing 
outwards. Here we must emphasize that orbits with initial 
conditions outside the interior region, or in other words out¬ 
side Li or L 2 , does not necessarily escape immediately from 
the galaxy. In Fig. we present a characteristic example 
of such an orbit with initial conditions: xq = —15.18 > tl, 
}^o = -^b = O 5 yo > 0 when C = 0.01. We observe that even 
though this orbit is initiated outside but very close to Li it 
does not escape right away from the system. On the other 
hand, it enters the interior galactic region and only after 32 
time units of chaotic motion it eventually escapes from the 
opposite channel. Thus it is evident, that the initial position 
itself does not furnish a sufficient condition for escape, since 
the escape criterion is a combination of the coordinates and 
the velocity of stars. 

The physical and the phase space are divided into the es¬ 
caping and non-escaping (trapped) space. Usually, the vast 
majority of the trapped space is occupied by initial con¬ 
ditions of regular orbits forming stability islands where a 
third integral is present. In many systems however, trapped 


chaotic orbits have also been observed. Therefore, we de¬ 
cided to distinguish between regular and chaotic trapped mo¬ 
tion. Over the years, several chaos indicators have been de¬ 
veloped in order to determine the character of orbits. In our 
case, we chose to use the Smaller ALingment Index (SALI) 
method. The SALI (Ml has been proved a very fast, reli¬ 
able and effective tool, which is defined as 

SALI(t) =min(d_,d+), (9) 

where J_ = ||wi(r) — W 2 ( 0 II and = ||wi(r) + W 2 ( 0 II are 
the alignments indices, while wi(r) and yv 2 {t), are two devi¬ 
ation vectors which initially point in two random directions. 
For distinguishing between ordered and chaotic motion, all 
we have to do is to compute the SALI along time interval 
tmax of numerical integration. In particular, we track simul¬ 
taneously the time-evolution of the main orbit itself as well 
as the two deviation vectors wi (t) and W 2(0 in order to com¬ 
pute the SALT 

The time-evolution of SALI strongly depends on the na¬ 
ture of the computed orbit since when an orbit is regular 
the SALI exhibits small fiuctuations around non zero val¬ 
ues, while on the other hand, in the case of chaotic orbits 
the SALI after a small transient period it tends exponen¬ 
tially to zero approaching the limit of the accuracy of the 
computer (10“^^). Therefore, the particular time-evolution 
of the SALI allow us to distinguish fast and safely between 
regular and chaotic motion. Nevertheless, we have to de¬ 
fine a specific numerical threshold value for determining the 
transition from order to chaos. After conducting extensive 
numerical experiments, integrating many sets of orbits, we 
conclude that a safe threshold value for the SALI is the value 
10“^. Thus, in order to decide whether an orbit is regular 
or chaotic, one may follow the usual method according to 
which we check after a certain and predefined time interval 
of numerical integration, if the value of SALI has become 
less than the established threshold value. Therefore, if SALI 
< 10“^ the orbit is chaotic, while if SALI > 10“^ the orbit 
is regular thus making the distinction between regular and 
chaotic motion clear and beyond any doubt. For the compu¬ 
tation of SALI we used the LP-VI code (24l, a fully oper¬ 
ational routine which efficiently computes a suite of many 
chaos indicators for dynamical systems in any number of 
dimensions. 

In our computations, we set 10"^ time units as a max¬ 
imum time of numerical integration. The vast majority of 
orbits (regular and chaotic) however, need considerable less 
time to find one of the two exits in the limiting curve and 
eventually escape from the system (obviously, the numer¬ 
ical integration is effectively ended when an orbit passes 
through one of the escape channels and escapes). Neverthe¬ 
less, we decided to use such a vast integration time just to 
be sure that all orbits have enough time in order to escape. 
Remember, that there are the so called “sticky orbits" which 
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Table 1 Type, initial conditions, value of the energy and integration 
time of the orbits shown in Eig.[^a-d). When the initial value of x or 
y is indicated as > 0 it means that is obtained from the Jacobi integral 

0 - 


Figure 

Type 

Vo 

yo 

Vo 

yb 

c 

^int 

3 

i 

1:1 loop 

-4.10 

0.00 

0 

>0 

0.01 

100 

3 

i 

1:1 loop 

-1.85 

0.00 

0 

>0 

0.01 

100 

3 


1:3 box 

0.00 

-1.25 

>0 

0 

0.10 

100 

3 


chaotic 

-1.00 

0.00 

>0 

0 

0.01 

164 

3 

1 

chaotic 

-2.40 

0.00 

>0 

0 

0.01 

129 


behave as regular ones during long periods of time. Here we 
should clarify, that orbits which do not escape after a numer¬ 
ical integration of 10 ^ time units ( 10 ^^ yr) are considered as 
non-escaping or trapped. In fact, orbits with escape periods 
equal to many Hubble times are completely irrelevant to our 
investigation since they lack physical meaning. 

4 Numerical results 

The main objective of our investigation is to determine which 
orbits escape and which remain trapped, distinguishing si¬ 
multaneously between regular and chaotic trapped motioij^ 
Furthermore, two additional properties of the orbits will be 
examined: (i) the channels through which the stars escape 
and (ii) the time-scale of the escapes (we shall also use the 
terms escape period or escape rates). In the present paper, 
we shall explore these dynamical quantities for various val¬ 
ues of the energy, always within the interval C G [0.001,0.1]. 

Our numerical calculations indicate that apart from the 
escaping orbits there is always a considerable amount of 
non-escaping orbits. In general terms, the majority of non¬ 
escaping regions corresponds to initial conditions of regular 
orbits, where a third integral of motion is present, restrict¬ 
ing their accessible phase space and therefore hinders their 
escape. However, there are also chaotic orbits which do not 
escape within the predefined interval of 10 ^ time units and 
remain trapped for vast periods until they eventually escape 
to infinity. At this point, it should be emphasized and clari¬ 
fied that these trapped chaotic orbits cannot be considered, 
by no means, neither as sticky orbits nor as super sticky or¬ 
bits with sticky periods larger than 10^ time units. Sticky 
orbits are those who behave regularly for long time periods 
before their true chaotic nature is fully revealed. In our case 
on the other hand, this type of orbits exhibit chaoticity very 
quickly as it takes no more than about 100 time units for the 

^ Generally, any dynamical method requires a sufficient time inter¬ 
val of numerical integration in order to distinguish safely between or¬ 
dered and chaotic motion. Therefore, if the escape rate of an orbit is 
very low or even worse if the orbit escapes directly from the system 
then, any chaos indicator (the SALI in our case) will fail to work prop¬ 
erly due to insufficient integration time. Hence, it is pointless to speak 
of regular or chaotic escaping orbits. 


SALI to cross the threshold value (SALI <C 10“^), thus iden¬ 
tifying beyond any doubt their chaotic character. Therefore, 
we decided to classify the initial conditions of orbits in both 
the physical and phase space into four main categories: (i) 
orbits that escape through Li, (ii) orbits that escape through 
L 2 , (iii) non-escaping regular orbits and (iv) trapped chaotic 
orbits. 

Additional numerical computations reveal that the non¬ 
escaping regular orbits are mainly 1:1 loop orbits for which 
a third integral applie^ while other types of secondary res¬ 
onant orbits are also present. In Fig.|^ we present two 1:1 
thin loop orbits, while in Fig.[^ a typical example of a sec¬ 
ondary 1:3 resonant quasi-periodic orbit is shown. The n : m 
notation we use for the regular orbits is according to 1231 
I112B . where the ratio of those integers corresponds to the 
ratio of the main frequencies of the orbit, where main fre¬ 
quency is the frequency of greatest amplitude in each coor¬ 
dinate. Main amplitudes, when having a rational ratio, define 
the resonances of an orbit. Finally in Figs.[^-d we observe 
two orbits escaping through Li (exit channel 1 ) and L 2 (exit 
channel 2), respectively. All regular orbits shown in Figs, 
[^-b were computed until /^ = 100 time units, while on the 
other hand, the escaping orbits presented in Figs.[^-d were 
calculated for 5 time units more than the corresponding es¬ 
cape period in order to visualize better the escape trail. In 
Table we provide the type, the exact initial conditions and 
the value of the energy for all the depicted orbits. 

A simple qualitative way for distinguishing between reg¬ 
ular and chaotic motion in Hamiltonian systems is by plot¬ 
ting the successive intersections of the orbits using a Poincare 
Surface of Section (PSS) ||49l. In particular, a PSS is a two- 
dimensional (2D) slice of the entire four-dimensional (4D) 
phase space. The following Fig. [^a-b) shows two PSSs at 
the critical Jacobi energy Cl. Fig.|^ corresponds to the phys¬ 
ical (v,y) space showing orbits crossing y = 0 with i > 0, 
while orbits crossing y = 0 with y > 0 at the phase (x,i) 
space are presented in Fig.|^. In both cases, the initial con¬ 
ditions of the orbits are integrated forwards in time plot¬ 
ting a dot at each crossing through the surface of section. 
The outermost black thick curve circumscribing each PSS is 
the limiting curve which in the physical space is defined as 
= E, while in the phase space is given by 

= o)=E. (10) 

It is seen, that the vast majority of the surfaces of section is 
covered by a unified and extended chaotic sea however, sta¬ 
bility regions with initial conditions corresponding to regu¬ 
lar orbits are also present. In these stability islands the quasi- 
periodic orbits admit a third integral of motion, also known 
as “adelphic integral" (see e.g., EZlX which hinders the test 
particles (stars) from escaping. Moreover, we observe that 

^ The total angular momentum is an approximately conserved quan¬ 
tity (integral of motion), even for orbits in non spherical potentials. 
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Fig. 3 Characteristic examples of the main types of orbits in our galaxy model, (a-upper left): two 1:1 loop orbits, (b-upper right): a 1:3 resonant 
quasi-periodic box orbit, (c-lower left): an orbit escaping through L \, (d-lower right): an orbit escaping through L 2 . More details are given in Table 

□ 




Fig. 4 Poincare Surfaces of Section at the critical Jacobi energy C^. (a-left): The physical (x,y) space where orbits crossing 3 ) = 0 with i: > 0 and 
(b-right): the phase (x,x) space where orbits crossing y = 0 with y > 0 . 
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some areas in the chaotic sea of the plane are less 
densely occupied than others. We remark that in the follow¬ 
ing grid exploration in the physical and phase space, we do 
not consider all initial conditions in the same chaotic domain 
as one the same orbit. 

4.1 The structure of physical (x,y) space 

Our exploration begins in the physical space and in the left 
column of Fig. 1^ we present the orbital structure of the (x,y) 
plane for three values of the energy, where the initial condi¬ 
tions of the orbits are classified into four categories by using 
different colors. Specifically, gray color corresponds to reg¬ 
ular non-escaping orbits, blue color corresponds to trapped 
chaotic orbits, green color corresponds to orbits escaping 
through channel 1, while the initial conditions of orbits es¬ 
caping through exit channel 2 are marked with red color. The 
vertical, dashed, magenta lines indicate the position of the 
two Lagrangian points Li and L 2 . For C = 0.001, that is an 
energy level just above the critical escape energy Cl, we see 
that the vast majority of initial conditions corresponds to es¬ 
caping orbits however, a stability island of non-escaping reg¬ 
ular orbits is present at the central upper part of the interior 
galactic region. We also see that the entire interior region 
is completely fractal which means that there is a highly de¬ 
pendence of the escape mechanism on the particular initial 
conditions of the orbits. In other words, a minor change in 
the initial conditions has as a result the star to escape through 
the opposite exit channel, which is of course, a classical in¬ 
dication of chaotic motion. With a much closer look at the 
(x,y) plane we can identify some additional tiny stability is¬ 
lands which are deeply buried in the escape domain. As we 
proceed to higher energy levels the fractal area reduces and 
several basins of escape start to emerge. By the term basin of 
escape, we refer to a local set of initial conditions that corre¬ 
sponds to a certain escape channel. Indeed for C = 0.01 we 
observe the presence of escape basins which mainly have 
the shape of thin elongated bands. Moreover, the extent of 
the central stability island seems to be unaffected by the in¬ 
crease in the orbital energy. With increasing energy the in¬ 
terior region in the physical spaces becomes less and less 
fractal and broad, well-defined basins of escape dominate 
when C = 0.1. The fractal regions on the other hand, are con¬ 
fined mainly near the boundaries between the escape basins 
or in the vicinity of the stability islands. Furthermore, it is 
seen that for C = 0.1 apart from the central stability island, 
there is also a set of three smaller islands of regular motion 
corresponding to 1:3 resonant orbits. It should be empha¬ 
sized that our classification shows that in all three energy 
levels trapped chaotic motion is almost negligible as the ini¬ 
tial conditions of such orbits appear only as lonely points 
around the boundaries of the stability islands. Here we must 
point out three interesting phenomena that take place with 



Fig. 6 Evolution of the percentages of escaping and non-escaping (reg¬ 
ular plus chaotic) orbits on the physical (x, y) space when varying the 
dimensionless energy parameter C. 


increasing energy: (i) the regions of forbidden motion are 
significantly reduced, (ii) the two escape channels become 
more and more wide and (iii) the amount of orbits with ini¬ 
tial conditions outside the interior region (especially outside 
L\) that escape from the opposite exit channel increases. 

In the right column of Fig. we show how the escape 
times 4sc of orbits are distributed on the physical (x,y) space. 
The escape time ^sc is defined as the time when a star passes 
one of the two vertical, dashed lines at (x,y) = (±rL,y) with 
velocity pointing outwards. Light reddish colors correspond 
to fast escaping orbits with short escape periods, dark blue/purpe 
colors indicate large escape rates, while white color denote 
both non-escaping regular and trapped chaotic orbits. The 
scale on the color bar is logarithmic. It is evident, that orbits 
with initial conditions close to the boundaries of the stabil¬ 
ity islands need significant amount of time in order to escape 
from the galaxy, while on the other hand, inside the basins 
of escape where there is no dependence on the initial con¬ 
ditions whatsoever, we measured the shortest escape rates 
of the orbits. We observe that for C = 0.001 the escape pe¬ 
riods of orbits with initial conditions in the fractal region 
are huge corresponding to tens of thousands of time units. 

As the value of the total orbital energy increases however, 
the escape times of orbits reduce significantly. In fact for 
C = 0.01 and C = 0.1 the basins of escape can be distin¬ 
guished in the grids of Fig.|^ being the regions with reddish 
colors indicating extremely fast escaping orbits. Our numer¬ 
ical calculations indicate that orbits with initial conditions 








Escape dynamics and fractal basin boundaries in Seyfert galaxies 


9 



X 

Fig. 5 Left column: Orbital structure of the physical (x,y) space. Top row: C = 0.001; Middle row: C = 0.01; Bottom row: C = 0.1. The green 
regions correspond to initial conditions of orbits where the stars escape through Li, red regions denote initial conditions where the stars escape 
through L 2 , gray areas represent stability islands of regular non-escaping orbits, initial conditions of trapped chaotic orbits are marked in blue, 
while the regions of forbidden motion are shown in yellow. The vertical, dashed, magenta lines indicate the position of the two Lagrangian points 
Li and L 2 . Right column: Distribution of the corresponding escape times fesc of the orbits on the physical space. The darker the color, the larger 
the escape time. The scale on the color-bar is logarithmic. Initial conditions of non-escaping regular orbits and trapped chaotic orbits are shown in 
white. 
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inside the basins have significantly small escape periods of 
less than 10 time units. 

The following Fig. [^presents the evolution of the per¬ 
centages of the different types of orbits on the physical (v,y) 
space when the dimensionless energy parameter C varies. 
At this point we must clarify that we decided to merge the 
percentages of non-escaping regular and trapped chaotic or¬ 
bits together because our computations indicate that always 
the rate of trapped chaotic orbits is extremely small (less 
than 0.1%) and therefore it does not contribute to the over¬ 
all orbital structure of the dynamical system. One may ob¬ 
serve, that when C = 0.001, that is just above the critical 
escape energy Cl, escaping orbits through Li and L 2 share 
about 96% of the available physical space, so the two exit 
channels can be considered equiprobable. As the value of 
the energy increases however, the percentages of escaping 
orbits through exit channels 1 and 2 start to diverge, espe¬ 
cially for C > 0.01. In particular, the amount of orbits that 
escape through Li exhibits a linear increase, while on the 
other hand, the rate of escaping orbits through the opposite 
exit channel displays a linear decrease. At the highest en¬ 
ergy level studied (C = 0.1) escaping orbits through Li is the 
most populated family occupying about 55% of the physical 
plane, while only about 40% of the same plane corresponds 
to initial conditions of orbits that escape through L 2 . Further¬ 
more, the percentage of non-escaping (regular plus chaotic) 
orbits is much less affected by the shifting of the energy dis¬ 
playing a rather constant value around 5%. Thus taking into 
account all the above-mentioned analysis we may conclude 
that at low energy levels where the fractility of the physi¬ 
cal space is maximum the stars do not show any particular 
preference regarding the escape channel, while on the con¬ 
trary, at high enough energy levels where basins of escape 
dominate it seems that exit channel 2 is more preferable. 

4.2 The structure of the phase {x^x) space 

We continue our investigation in the phase (v,i) space and 
we follow the same numerical approach as discussed previ¬ 
ously. In Fig. 1^ we depict the orbital structure of the (v,i) 
plane for three values of the energy, using different colors 
in order to distinguish between the four main types of orbits 
(non-escaping regular; trapped chaotic; escaping through Li 
and escaping through L 2 ). Just above the escape energy, that 
is when C = 0.001, we see that the structure of the phase 
plane is very similar to that of the corresponding physical 
plane shown in Fig. It is observed that about 90% of 
the (v,i) plane is occupied by initial conditions of escaping 
orbits. Moreover, the vast escape domain is highly fractal, 
while basins of escape, if any, are negligible. In all studied 
cases the areas of regular motion correspond mainly to retro¬ 
grade orbits (i.e., when a star revolves around the galaxy in 
the opposite sense with respect to the motion of the galaxy 


100 



X 


Fig. 8 Reconstruction of the phase (x,i) plane for C = 0.1 where the 
color scale of the escape regions is given as a function of the number 
of intersections with the = 0 axis upwards ( 3 ) > 0). The color code is 
as follows: 0 intersections (red); 1 intersection (blue); 2 intersections 
(magenta); 3-10 intersections (orange); >10 intersections (green). The 
gray regions represent stability islands of non-escaping orbits. 

itself), while there are also some smaller stability islands 
of prograde orbits at high energy levels (i.e., for C = 0.1). 
The area on the phase plane covered by escape basins grows 
drastically with increasing energy and at high enough energy 
levels the fractal regions are confined mainly at the bound¬ 
aries of the stability islands of non-escaping orbits. In par¬ 
ticular, the escape basins first appear at mediocre values of 
energy (C = 0.01) inside the fractal, extended escape region 
having the shape of thin elongated bands however, as we 
proceed to higher energy levels they grow in size dominating 
in the phase space as well-formed broad domains. It is worth 
noticing the initial conditions of trapped chaotic orbits at the 
central region of the stability islands for C = 0.1 which indi¬ 
cates the presence of a thin chaotic layer (separatri}0. The 
Coriolis forces dictated by the rotation of the compact active 
nucleus of the Seyfert galaxy makes the phase planes to be 
asymmetric with respect to the i-axis and this phenomenon 
is usually known as “Coriolis asymmetry" (e.g., 1521 ). 

The distribution of the corresponding escape times 4sc 
of orbits on the phase space as a function of the energy pa¬ 
rameter is shown in the right column Fig. |7] One may ob¬ 
serve that the results are very similar to those presented ear¬ 
lier in Fig. where we found that orbits with initial con¬ 
ditions inside the basins of escape have the smallest escape 


^ The separatrix prevents chaotic orbits from escaping in the same 
manner that an adelphic integral hinders quasi-periodic orbits from es¬ 
caping. 
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Fig. 7 Left column: Orbital structure of the phase (x,i) space. Top row: C = 0.001; Middle row: C = 0.01; Bottom row: C = 0.1. Right column: 
Distribution of the corresponding escape times ^esc of the orbits on the phase space. The color codes are the same as in Fig.|^ 


rates, while on the other hand, the longest escape times cor¬ 
respond to orbits with initial conditions either in the fractal 
regions of the plots, or near the boundaries of the islands 


of regular motion. Another interesting way of measuring the 
escape period of an orbit is by counting how many intersec¬ 
tion the orbit has with the axis y = 0 before it escapes. The 
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Fig. 9 Evolution of the percentages of escaping and non-escaping (reg¬ 
ular plus chaotic) orbits on the phase (x,i:) space when varying the 
dimensionless energy parameter C. 
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regions of the phase space for C = 0.1 in Fig. are colored 
according to the number of intersections with the axis y = 0 
upwards (y > 0). Our computations reveal that orbits that 
escape directly without any intersection with the y = 0 axis 
(0 intersections) or performing between 3 and 10 intersec¬ 
tions or more than 10 intersections share about three fourths 
of the escape region. Orbits that cross one time the y = 0 
axis before escape correspond to about 15% of the escape 
realm, while only about 10% of the escaping orbits perform 
two intersections before leaving the system passing through 
one of the Lagrangian points. It is evident, that orbits with 
initial conditions located near the boundaries of the stabil¬ 
ity islands perform numerous intersections (> 10) with the 
y = 0 axis before they eventually escape to infinity. On the 
contrary, orbits with initial conditions in the escape basins 
need only a couple of intersection until they escape. 

The evolution of the percentages of the different types 
of orbits on the phase (x^x) space as a function of the di¬ 
mensionless energy parameter C is presented in Fig.[^ Once 
more, we have to point out that the percentages of non¬ 
escaping ordered and trapped chaotic orbits are merged to¬ 
gether in a single trend line because the percentage of trapped 
chaotic orbits is extremely small (less than 0.5%) through¬ 
out. It is seen, that at the lowest energy level studied (C = 
0.001) escaping orbits through Li and L 2 share about 90% 
of the phase plane, while only 10% of the same plane cor¬ 
responds to initial conditions of non-escaping (regular plus 
chaotic) orbits. Furthermore, we observe that for 0.001 < 
C < 0.01 the rates of all types of orbits remain practically 


unperturbed by the shifting on the value of the orbital en¬ 
ergy. For larger values of the energy (C > 0.01) however, the 
percentages of escaping orbits exhibit a small divergence, 
while the portion of non-escaping orbits displays a minor 
increase and for C = 0.1 they reach about 12.5%. Specifi¬ 
cally, the rate of orbits escaping through Li slightly reduces, 
while that of orbits escaping through L 2 continues the mono¬ 
tone behavior. It should be emphasized that the divergence 
of the percentages of escaping orbits observed in the phase 
space is much smaller than that found to exist in the physi¬ 
cal space. Therefore, one may reasonably deduce that in the 
phase space the increase in the value of the energy does not 
infiuence significantly the orbital content of the dynamical 
system, thus leaving the percentages of the orbits almost the 
same throughout. In addition, since that the rates of escap¬ 
ing orbits remain unperturbed within the energy range, we 
may say that the two exit channels can be considered almost 
equiprobable in the phase space. 

4.3 An overview analysis 

The color-coded grids in the physical (v,y) as well as in 
the phase (v, x) space provide sufficient information on the 
phase space mixing however, for only a fixed value of the 
Jacobi constant. Henon back in the late 60s Ell, introduced 
a new type of plane which can provide information not only 
about stability and chaotic regions but also about areas of 
trapped and escaping orbits using the section y = i = 0,y>0 
(see also CD- In other words, all the orbits of the stars of 
the galaxy are launched from the x-axis with x = xq, par¬ 
allel to the y-axis (y = 0). Consequently, in contrast to the 
previously discussed types of planes, only orbits with peri- 
centers on the x-axis are included and therefore, the value 
of the dimensionless energy parameter C can be used as an 
ordinate. In this way, we can monitor how the energy infiu- 
ences the overall orbital structure of our dynamical system 
using a continuous spectrum of energy values rather than 
few discrete energy levels. In Fig. p^ we present the orbital 
structure of the (x,C)-plane when C G [0.001,1], while in 
Fig. \T0j^ the distribution of the corresponding escape times 
of orbits is depicted. 

We observe that for relatively low energy values (0.001 < 
C < 0.01) the interior galactic region is highly fractal, while 
basins of escape are located only at the outer parts of the 
(x,C)-plane, that is outside Li and L 2 (or in other words in 
the exterior galactic region). Furthermore, we can identify 
the presence of two small stability islands of prograde (xq > 
0) non-escaping regular orbits which were invisible in the 
general plots shown earlier in Figs.[^and[7] For larger values 
of energy C > 0.01 however, the structure of the (x, C)-plane 
changes drastically and the most important differences are 
the following: (i) several basins of escape are formed inside 
the fractal escape region, (ii) apart from the main stability 
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Fig. 10 Orbital structure of the (a-upper left): (x,C)-plane; and (b-upper right): (y,C)-plane; (c-lower left and d-lower right): the distribution of 
the corresponding escape times of the orbits. The color codes are the exactly same as in Fig.[^ 


island a second smaller one emerges near Li. It should be 
pointed out that in the blow-ups of the diagram several ad¬ 
ditional extremely tiny islands of stability have been identi- 
fiecj^ (hi) at high enough energy levels (C > 0.1) we see that 
a thin chaotic layer composed of initial conditions of trapped 
chaotic orbits appear at the central region of the main stabil¬ 
ity island and (iv) the extent of the main island of regular 
motion grows rapidly and for C > 1 the position of initial 
conditions of non-escaping regular orbits exceed the interior 
region {xq > vl). The following Tableshows the percent¬ 
ages of the four main types of orbits in the color-coded grids 
shown in Fig. 10 a-b). 


In order to obtain a more complete view on how the total 
orbital energy influences the nature of orbits in our galaxy 
model, we follow a similar numerical approach to that ex- 


^ An infinite number of regions of (stable) quasi-periodic (or small 
scale chaotic) motion is expected from classical chaos theory. 


Table 2 Percentages of the four main types of orbits in the color-coded 
grids shown in Figs.p^a-b). 


Figure 

Trapped 

Non-escaping 

Channel 1 

Channel 2 

10 

i 

0.34 

19.05 

31.94 

48.67 

10 


1.46 

28.72 

26.53 

43.29 


plained before but now all orbits of stars are initiated from 
the vertical y-axis with y = yo- In particular, this time we use 
the section x = y = 0 , i > 0 , launching orbits parallel to the 
v-axis. This allow us to construct again a 2D plane in which 
the y coordinate of orbits is the abscissa, while the logarith¬ 
mic value of the energy logiQ(C) is the ordinate. The orbital 
structure of the (y,C)-plane when C G [0.001,1] is shown 
in Fig.p^. The vertical, dashed, magenta lines indicate the 
position of the Lagrangian points Lq and L 5 . The black solid 
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line is the limiting curve which distinguishes between re¬ 
gions of allowed and forbidden motion and is defined as 


/L(j,C) = 0eff(O,^)=^. 


( 11 ) 


A very complicated orbital structure is reveled in the (y,C)- 
plane which however, in general terms, is very similar with 
that of the (v, C)-plane. One may observe that for E > —374.55 
the forbidden regions of motion around L 4 and L 5 completely 
disappear. In fact the value E{L/[) = E{Ls) = —374.55 is 
another critical energy level. It is seen that above that crit¬ 
ical value the left part (yo < 0 ) of the plane is occupied 
almost completely with initial conditions of orbits that es¬ 
cape through L 2 , while the right part of the same plane on 
the other hand, contains mostly initial conditions of non¬ 
escaping regular orbits. 


It is evident from the results presented in Figs. 10 c-d) 
that the escape times of the orbits are strongly correlated 
to the escape basins. In addition, one may conclude that the 
smallest escape periods correspond to orbits with initial con¬ 
ditions inside the escape basins, while orbits initiated in the 
fractal regions of the planes or near the boundaries of sta¬ 
bility islands have the highest escape rates. In both types of 
planes the escape times of orbits are significantly reduced 
with increasing energy. Combining all the numerical out¬ 
comes presented in Figs. [T^a-d) we may say that the key 
factor that determines and controls the escape times of the 
orbits is the value of the orbital energy (the higher the en¬ 
ergy level the shorter the escape rates), while the fractality 
of the basin boundaries varies strongly both as a function of 
the energy and of the spatial variable. 

The evolution of the average value of the escape time 
< 4 sc > of orbits as a function of the dimensionless energy 
parameter is given in Fig.pT^. It is seen, that for low values 
of energy the average escape time of orbits reduces rapidly, 
while for C > 0.01, it seems to saturate around 80 time units. 
We feel it is important to justify this behaviour of the escape 
time. As previously explained, a Lyapunov orbit acts like a 
barrier, reaches the ZVC on both sides of the opening and 
returns along the same path thus, connecting two opposite 
branches of the ZVC. Therefore, we could exploit this fact 
in order to quantify the escape process. For this purpose, we 
define w as the distance between the two points at which 
the Lyapunov orbit touches the two opposite branches of the 
ZVC. Using w we can measure, in a way, the width of the 
escape channels. Fig. \n\ ) depicts the evolution of w as a 
function of the energy parameter C. We observe an almost 
exponential growth of w with increasing energy. This means 
that as the value of the Jacobi constant increases the escape 
channels (which are of course symmetrical) become more 
and more wide and therefore, orbits need less and less time 
in order to find one of the two openings in the open equipo- 
tential surface and eventually escape from the system. This 
geometrical feature explains why for low values of energy 
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Fig. 12 Visualization of the spiral path of two orbits after escaping 
from Li (green) and L 2 (red) when C = 0.01. The blue dots indicate the 
position of Li and L 2 , while the yellow banana-shaped isolines define 
the forbidden regions of motion around L 4 and L 5 . 


orbits consume large time periods wandering inside the open 
equipotential surface until they eventually locate one of the 
two exits and escape to infinity. Therefore, we may conclude 
that it is the geometrical shape of the exit channels that jus¬ 
tifies the decreasing pattern of ^esc shown in Fig.pT^. 


5 Formation of spiral arms 

It is well known that when stars escape from star clusters 
through the Lagrangian points form complicated structures 
known as tidal tails or tidal arms (e.g., 01381 EH E21IM1)- 
In the same vein, in the case of barred spiral galaxies we 
observe the formation of spiral arms (e.g., 144114^17411751 
[851 ESI |90l). Usually a star cluster rotates around its par¬ 
ent galaxy in a circular orbit thus, the tidal tails follow the 
curvature defined by the circular orbit and they are nearly 
straights for orbits moving in large galactocentric distances 
(see e.g., ED). The spiral arms of barred galaxies on the 
other hand, are formed by the non-axisymmetric perturba¬ 
tion of the bar. In particular, in barred spiral galaxies with 
prominent spiral-like shape the two arms start from the ends 
of the bar and then they wind up around the banana-shaped 
forbidden regions which enclose the Lagrangian points L 4 
and L 5 . 

Almost all Seyfert galaxies are spiral galaxies and have 
been among the most intensively studied objects in astron¬ 
omy (e.g., 1891). Therefore, it would be very challenging to 
investigate if our simple gravitational model has the ability 
to realistically simulate the formation of spiral arms. As a 
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Fig. 11 Evolution of (a-left): the average escape time of orbits < ^esc > and (b-right): the width w of the escape channels as a function of the orbital 
energy C. 


first step, we integrated the two orbits shown in Fig. [^c-d) 
for 7 time units after crossing the Lagrangian points and in 
Fig. [T^ we present the corresponding paths of the orbits that 
escape from opposite exits. It is seen, that both orbits fol¬ 
low indeed a spiral path around the yellow banana-shaped 
isolines of the effective potential in which the Lagrangian 
points L 4 and L 5 are enclosed. 


Taking into account that our initial tests are positive re¬ 
garding the formation of the spiral structure, we reconfig¬ 
ured our integration routine for the computation of the basins 
of escape to yield output of all orbit trajectories for a three- 
dimensional grid of size x x = 100 x 100 x 100 for 
C = 0.01. This dense grid of initial conditions is centered 
at the origin of the phase space and we allow for all orbits 
both sings of y. Our purpose is to monitor the time-evolution 
of the path of the orbits and simulate the spiral arm forma¬ 
tion. At t = 0 all orbits are regularly distributed within the 
Lagrangian radius inside the interior galactic region. In Fig. 
p^we present four time snapshots of the position of the 10 ^ 
stars in the physical (v,y) plane. We observe that aXt = 1.25 
time units, which corresponds to 125 Myr the majority of 
stars are still inside the interior galactic region however, a 
small portion of stars has already escaped passing through 
either Li or L 2 . At ^ = 2.5 time units (250 Myr) we have the 
first indication of the formation of spiral arms, where with 
green and red color we depict stars that escaped through exit 
channels 1 and 2 , respectively, while stars that remain inside 
the Lagrangian radius are shown in gray. As time goes by we 
see that the two symmetrical spiral arms grow in size and at 
t = 5 time units (500 Myr) the Seyfert galaxy has obtained 
its complete spiral structure, while the interior galactic re¬ 
gion is uniformly filled. At this point we have to point out 


that the standard integration code, which was used to obtain 
the results presented in Figs. and [7] plots a point at ev¬ 
ery step of the numerical integration. In Fig.[^a-d) on the 
other hand, the density of points along one star orbit is taken 
to be proportional to the velocity of the star. Being more 
specific, a point is plotted (showing the position of a star), if 
an integer counter variable which is increased by one at ev¬ 
ery integration step, exceeds the velocity of the star. Follow¬ 
ing this technique we can simulate, in a way, a real A-body 
simulation of the spiral evolution of the galaxy, where the 
density of stars will be highest where the corresponding ve¬ 
locity is lowest. Therefore, we proved that for C = 0.01 the 
scenario of evolution of spiral arms in our galaxy model is 
indeed viable. Additional numerical simulations reveal that 
this is also true for other (lower or higher) values of energy. 
In fact, the particular energy level mainly controls how tight 
the spiral arms wind up around the forbidden regions of mo¬ 
tion. 

In section]^ we explained that the values of the mass Mg 
and the scale length c entering Eq. Q were chosen having 
in mind the Seyfert galaxy NGC 1566 with total mass equal 
to 1.2 X10^^ M 0 1^ . Here we have to point out that NGC 
1566 is not just any Seyfert galaxy; it is the second brightest 
Seyfert galaxy known. It is also the brightest and most dom¬ 
inant member of the Dorado Group, a loose concentration of 
galaxies located approximately 40 million light-years away 
that together comprise one of the richest galaxy groups of 
the southern hemisphere. NGC 1566 is an intermediate spi¬ 
ral galaxy of type SAB(rs)bc, meaning that while it does not 
have a well defined bar-shaped region of stars at its central 
region - like normal barred spirals - it is not quite an un¬ 
barred spiral either. The small but extremely bright nucleus 
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Fig. 13 Simulation time snapshots of the position of 10^ stars in the physical (x,y) plane initiated {t = 0) within the Lagrangian radius ri, for 
C = 0.01. We see that as time evolves two symmetrical spiral arms are formed. The green arm contains stars that escaped through Li, while the 
red arm contains stars that escaped through L 2 . The bounded (trapped or non-escaping) stars inside the interior galactic region are shown in gray. 


of NGC 1566 is clearly visible in Fig. II105L a telltale 
sign of its membership of the Seyfert class of galaxies. The 
centres of such galaxies are very active and luminous, emit¬ 
ting strong bursts of radiation and potentially harbouring 
super-massive black holes that are many millions of times 
the mass of the Sun. This image highlights the beauty and 
awe-inspiring nature of this unique galaxy group, with NGC 
1566 glittering and glowing, its bright nucleus framed by 
swirling and symmetrical lavender arms. In Fig. we have 
superposed the time snapshot ait = 5 time units (500 Myr) 
of our simulation above the real image of the galaxy. We ob¬ 
serve that the central interior region adequately matches the 
active galactic core, while the two spiral arms obtained from 
the numerical simulation fit the corresponding real ones. 
Thus, we may conclude that our simple one-component 
gravitational effective potential Q can, in a way, model not 


only the escape dynamics but also the spiral formation and 
evolution of Seyfert galaxies. 

Undoubtedly, a Seyfert galaxy is a very complex stel¬ 
lar entity and therefore, we need to assume some neces¬ 
sary simplifications and assumptions in order to be able to 
study mathematically the orbital behavior of such a com¬ 
plicated stellar system. For this purpose, our gravitational 
model is intentionally simple and contrived, in order to give 
us the ability to study all the different aspects regarding the 
kinematics and dynamics of the system. Nevertheless, con¬ 
trived models can surely provide an insight into more real¬ 
istic stellar systems, which unfortunately are very difficult 
to be explored, if we take into account all the astrophysi- 
cal aspects (i.e., gas, mergers, close encounters etc). Self- 
consistent models on the other hand, are mainly used when 
conducting W-body simulations. However, this application 
is entirely out of the scope of the present paper. Once more. 
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Fig. 14 A real image of the Seyfert galaxy NGC 1566 where we have 
superposed the simulation output of our gravitational model at r = 500 
Myr shown in Fig.pjjl. The gravitational potential studied in this paper 
is overall potential of the galaxy. 

we have to point out that the simplicity of our model is nec¬ 
essary, otherwise it would be extremely difficult, or even im¬ 
possible, to apply the extensive and detailed numerical cal¬ 
culations presented in this study. Our single-component po¬ 
tential neglect many features which must exist in real Seyfert 
galaxies such as colors/ages of stars populations, deviations 
from simple symmetries and time-independent perturbations. 
However, it is still significant that the conclusions derived 
from this simple potential all agree, both qualitatively and 
semi-quantitatively. This suggests that our model is robust, 
depending only on such topological properties of the phase 
space. 

6 Discussion and conclusions 

The aim of this work was to shed some light to the trapped or 
escaping nature of orbits in a dynamical gravitational model 
describing a Seyfert galaxy. In order to describe the proper¬ 
ties of the Seyfert galaxy as a whole we used an asymmetric 
Plummer-type potential. At this point we should emphasize 
that a Seyfert galaxy is composed of a disk, bulge a cen¬ 
tral nucleus and a bar. However, our dynamical system has 
two degrees of freedom and for this type of Hamiltonians 
we usually adopt single component potentials to model the 
properties of barred galaxies (e.g., lfT^fT^[T^I^[30l[36ll ) 
or galaxies in general (e.g., Il9l [^[^l^l63l[64l[73lll06ll ). 
This approach allow us to directly determine the effects of 
the single component potential to the the different types (res¬ 
onances) of orbits, the amount of chaos, the escape proper¬ 
ties and other dynamical aspects. 


The Jacobi integral of motion has a critical threshold 
value (or escape energy) which determines if the escape pro¬ 
cess is possible or not. In particular, for energies smaller than 
the threshold value, the equipotential surface is closed and 
it is certainly true that escape is impossible. For energy lev¬ 
els larger than the escape energy however, the equipotential 
surface opens and two exit channels appear through which 
the particles can escape to infinity. Here it should be em¬ 
phasized, that if a test particle does has energy larger than 
the escape value, there is no guarantee that the star will cer¬ 
tainly escape from the system and even if escape does occur, 
the time required for an orbit to transit through an exit and 
hence escape to infinity may be vary long compared with the 
natural escaping time. We managed to distinguish between 
ordered/chaotic and trapped/escaping orbits and we also lo¬ 
cated the basins of escape leading to different exit channels, 
finding correlations with the corresponding escape times of 
the orbits. Our extensive and thorough numerical investiga¬ 
tion strongly suggests, that the overall escape mechanism 
is a very complicated procedure and very dependent on the 
value of the Jacobi integral. 

Since a distribution function of the model was not avail¬ 
able so as to use it for extracting the different samples of 
orbits, we had to follow an alternative path. We defined for 
several values of the Jacobi integral, dense grids of initial 
conditions (vo,yo) and (vo,vb) regularly distributed in the 
area allowed by the energy level on the physical and phase 
space, respectively and then we identified regions of order/chaos 
and bound/escape. In each type of grid the step separation 
of the initial conditions along the axes (or in other words the 
density of the grid) was controlled in such a way that there 
are always at least 10^ orbits to be integrated and classified. 

For the numerical integration of the orbits in each type of 
grid, we needed about between 1 hour and 5 days of CPU 
time on a Pentium Dual-Core 2.2 GHz PC, depending on 
the escape rates of orbits in each case. For each initial con¬ 
dition, the maximum time of the numerical integration was 
set to be equal to 10"^ time units however, when a particle 
escaped the numerical integration was effectively ended and 
proceeded to the next available initial condition. 

The present article provides quantitative information re¬ 
garding the escape dynamics in spiral Seyfert galaxies. The 
main numerical results of our research can be summarized 
as follows: 

1. In all examined cases, areas of bounded motion and re¬ 
gions of initial conditions leading to escape in a given di¬ 
rection (basins of escape), were found to exist in both the 
physical and the phase space. The several escape basins 
are very intricately interwoven and they appear either as 
well-defined broad regions or thin elongated bands. Re¬ 
gions of bounded orbits first and foremost correspond to 
stability islands of regular orbits where a third integral 
of motion is present. 
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2. A strong correlation between the extent of the basins of 
escape and the value of the Jacobi integral was found to 
exist. Indeed, for low energy levels the structure of both 
the physical and the phase space exhibits a large degree 
of fractalization and therefore the majority of orbits es¬ 
cape choosing randomly exit channels. As the value of 
the energy increases however, the structure becomes less 
and less fractal and several basins of escape emerge. The 
extent of these basins of escape is more prominent at 
high energy levels. 

3. In both the physical and the phase space we identified a 
small portion of trapped chaotic orbits which do not es¬ 
cape within the predefined maximum time of numerical 
integration. The initial conditions of these orbits are lo¬ 
cated either in thin chaotic layers (separatrix) inside the 
regular areas of motion, or near the boundaries of sta¬ 
bility islands. Moreover, these orbits reveal their chaotic 
nature very quickly, so they cannot be considered as sticky 
orbits. 

4. We observed, that in many cases the escape process is 
highly sensitive dependent on the initial conditions, which 
means that a minor change in the initial conditions of an 
orbit lead the test particle to escape through the other exit 
channel. These regions are the exact opposite of the es¬ 
cape basins, are completely intertwined with respect to 
each other (fractal structure) and are mainly located in 
the vicinity of stability islands. This sensitivity towards 
slight changes in the initial conditions in the fractal re¬ 
gions implies, that it is impossible to predict through 
which exit the particle will escape. 

5. Our calculations revealed, that the escape times of orbits 
are directly linked to the basins of escape. In particular, 
inside the basins of escape as well as relatively away 
from the fractal domains, the shortest escape rates of the 
orbits had been measured. On the other hand, the longest 
escape periods correspond to initial conditions of orbits 
either near the boundaries between the escape basins or 
in the vicinity of the stability islands. 

6. It was detected, that for energy levels slightly above the 
escape energy the majority of the escaping orbits have 
considerable long escape rates (or escape periods), while 
as we proceed to higher energies the proportion of fast 
escaping orbits increases significantly. This behavior was 
justified, by taken into account a geometrical feature of 
the equipotential surface. Specifically, for low values of 
energy the width of the exit channels is too small there¬ 
fore, stellar orbits consume large time periods wandering 
inside the open equipotential surface until they eventu¬ 
ally locate one of the two exits and escape to infinity. 
Moreover, it was proved that the width of the escape 
channels grows rapidly with increasing energy, while at 
the same time the average escape time decreases signifi¬ 
cantly. 


7. We presented evidence that escaping stars with values of 
energy higher than the critical Jacobi value form spiral 
arms outside the interior galactic region defined by the 
Lagrangian radius. This procedure has a striking simi¬ 
larity to the formation of tidal tails in star clusters. We 
also conducted numerical simulations proving that our 
gravitational effective potential can model the forma¬ 
tion as well as the evolution of spiral arms in real Seyfert 
galaxies. 

Judging by the detailed and novel outcomes we may say 
that our task has been successfully completed. We hope that 
the present numerical analysis and the corresponding results 
to be useful in the field of escape dynamics in Active Galac¬ 
tic Nuclei (AGN). The outcomes as well as the conclusions 
of the present research are considered, as an initial effort 
and also as a promising step in the task of understanding the 
escape mechanism of orbits in AGNs. Taking into account 
that our results are encouraging, it is in our future plans to 
modify properly our galactic model in order to expand our 
investigation into three dimensions and explore the entire 
six-dimensional phase space. In addition, A-body simula¬ 
tions may elaborate the formation as well as the evolution of 
spiral structure in Seyfert galaxies. 
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